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Recently Y. N. showed that the nonequilibrium critical relaxation of the 2D Ising model from 
the perfectly-ordered state in the Wolff algorithm is described by the stretched-exponential decay, 
and found a universal scaling scheme to connect nonequilibrium and equilibrium behaviors. In the 
present study we extend these findings to vector spin models, and the 3D Heisenberg model could be 
a typical example. In order to evaluate the critical temperature and critical exponents precisely with 
the above scaling scheme, we calculate the nonequilibrium ordering from the perfectly-disordered 
state in the Swendsen-Wamg algorithm, and find that the critical ordering process is described by 
the stretched-exponential growth with the comparable exponent to that of the 3D XY model. The 
critical exponents evaluated in the present study are consistent with those in previous studies. 
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I. INTRODUCTION 

The cluster algorithm was first proposed [l[ in the 
Potts model [2!], which is a generalized version of the Ising 
model and characterized by short-range interactions be¬ 
tween discrete local variables. The partition function of 
this model can be represented with respect to percola¬ 
tion clusters [3j, and the cluster algorithm is formulated 
on the basis of this representation. In this algorithm 
temperature is included in the probability for connecting 
local variables to construct the clusters, and the flipping 
rate of clusters does not depend on temperature. Then, 
this algorithm has been expected to reduce the critical 
slowing down in the vicinity of the critical temperature 
T c , which is characterized by the dynamical critical ex¬ 
ponent 2 defined by divergence of the relaxation time of 
autocorrelation functions as t(T) ~ (T — T c )~ z . 

In finite systems the relaxation time is finite even at 
T c , and the exponent z is estimated from its size depen¬ 
dence as t(T c , L) ~ L z . Large reduction of z from that in 
the local-update algorithms was exhibited even in the pi¬ 
oneering articles of the cluster algorithms 0,1- In some 
numerical calculations it was pointed out that the weak 
power-law and logarithmic size dependences of t{L) were 
difficult to distinguish il- Such logarithmic size de¬ 
pendence stands for z = 0, namely the power-law dynam¬ 
ical critical behaviors may be questionable in the cluster 
algorithms. Then, much larger-scale calculations were 
performed by Giindiig et al. {fj and Du et al. Q with the 
nonequilibrium relaxation (NER), and they found that 
z = 0 widely holds in the cluster algorithms. 

Recently one of the present authors (Y. N.) found [9j] 
that the explicit form of NER of the magnetization in 
the 2D Ising model with the Wolff algorithm from the 
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perfectly-ordered state is described by the stretched- 
exponential time dependence, and proposed a univer¬ 
sal scaling scheme in which the whole relaxation process 
from early-time nonequilibrium behaviors to equilibrium 
ones for various system sizes is located on a single curve. 
In the present article we generalize these findings to vec¬ 
tor spin models, and choose the 3D Heisenberg model as 
a typical example, since it shows the second-order phase 
transition and the critical behaviors have already been 
studied intensively, similarly to the Ising models. 


Outline of the present article is as follows: In Sec¬ 
tion II, formulation based on the “embedded-Ising-spin” 
algorithm is briefly reviewed, and the reason why our 
calculations is based on the ordering from the perfectly- 
disordered state in the Swendsen-Wang (SW) algorithm 
is explained. In Section III, explicit scaling forms of var¬ 
ious physical quantities with size-dependent factors are 
exhibited. The critical temperature T c is evaluated from 
the ordering process of magnetization together with the 
critical exponent / 3 /u, and other exponents 7/1/ and a/i/ 
are estimated from the scaling behaviors of the magnetic 
susceptibility and specific heat. Then, the exponent v 
can be obtained from the scaling relation. However, eval¬ 
uation of ol/v from the specific heat is difficult in 3D vec¬ 
tor spin models, and v is also estimated from the scaling 
behavior of the temperature derivative of magnetization. 
In section IV, numerical results of the present articles are 
discussed. Estimates of the critical exponents /3, 7 and v 
are compared with other numerical studies, and promis¬ 
ing future tasks along the present framework are men¬ 
tioned. In section V, the above descriptions are summa¬ 
rized. Evaluation of T c from single system and nonequi¬ 
librium relaxation from the perfectly-ordered state in the 
SW and Wolff algorithms are treated in Appendices. 
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II. FORMULATION 

In the present article, Monte Carlo simulations of the 
ferromagnetic 3D Heisenberg model on a cubic lattice, 


When similar calculations are started from the perfectly- 
disordered state, the power-law ordering is observed in 
the Wolff algorithm [§|. On the other hand, such order¬ 
ing is described by the stretched-exponential one, 


n = -J J2 &-SJ (i) 

(ij')£n.n. 

with summation over all the nearest-neighbor bonds, are 
performed with the cluster algorithm. Wolff Q showed 
that the cluster update of vector spins such as Heisenberg 
one is possible by constructing spin clusters with respect 
to the Ising element projected onto a randomly-chosen 
direction r, |r| = 1 in each Monte Carlo step. That is, 
two nearest-neighbor spins S) and Sj are connected with 
probability p = 1 — exp[—2/3(r ■ Si)(f • 5j)] with /? = J/T, 
when the projected Ising elements along r, namely ( r-Si)r 
and (r • Sj)r, are aligned along the same direction. 

Wolff proposed 0) the so-called Wolff algorithm in 
the same article, in which a single spin cluster gener¬ 
ated from a randomly-chosen spin is always flipped. On 
the other hand, even in Wollff’s “embedded-Ising-spin” 
scheme, the SW update [l] is also possible, in which all 
spins are swept in every step and clusters with the flip¬ 
ping rate 50% are generated in the whole system. Al¬ 
though both approaches result in the same equilibrium, 
dynamical process of them might be different. Actu¬ 
ally, even in the Ising model, dynamics of these two al¬ 
gorithms are different in the ordering process from the 
perfectly-disordered state [§]. While physical quantities 
show stretched-exponential behaviors in the SW algo¬ 
rithm, they show power-law behaviors in the Wolff al¬ 
gorithm. 

While the cluster to be grown to system-size one is 
swept in each step in the SW algorithm by definition, 
such an event is quite rare in the Wolff algorithm. Since 
calculations from the perfectly-disordered state is re¬ 
quired for the NER analysis of physical quantities with 
respect to fluctuations [10, IHl ■ we utilize the SW algo¬ 
rithm in the present article. Furthermore, dynamical pro¬ 
cess in these two cases are different even in calculations 
from the perfectly-ordered state in the 3D Heisenberg 
model, which will be explained in Appendix B EH ■ 


III. NUMERICAL RESULTS 

A. Overview of calculations 

Recently one of the present authors (Y. N.) found @ 
that nonequilibrium critical relaxation of the absolute 
value of magnetization in the 2D Ising model is described 
by the stretched-exponential decay, 

<|m(f)|) ~ exp l-{t/T' m y) , 0 < cr < 1, (2) 

in the Wolff algorithm from the perfectly-ordered state. 
This scaling form also holds in the SW algorithm. 


(|m(t, L)\) ~ L d/ 2 exp[+(t/r m ) CT ], 0 < a < 1, (3) 

in the SW algorithm EH- Here early-time size depen¬ 
dence of the magnetization is also taken into account, 
and this form is derived from normalized random-walk 
growth of clusters. In the present article we also analyze 
other physical quantities on the basis of similar scaling 
forms. In order to investigate such diverging quantities 
at T c , calculations from the perfectly-disordered state are 
required in the NER method. Since the power-law order¬ 
ing in the Wolff algorithm has no essential difference from 
that in the local-update algorithm, here we concentrate 
on calculations based on the SW algorithm. 

The purpose to treat some physical quantities are to 
evaluate critical exponents. Although they do not appear 
explicitly in Eq. ©> the magnetization at the critical 
temperature T c converges to the critical one scaled as 
m c (L) ~ L~P/ V in equilibrium. Then, multiplying L^l v 
for both sides of Eq. © and taking logarithm, we have 

log((\m(t,L)\)L^) ~ (t/T m y -logL d / 2 -^. (4) 


When lhs of Eq. © is plotted versus rhs of it for 
several system sizes, short-time (in the region Eq. © 
holds) and long-time (in the vicinity of equilibrium) 
data are scaled by definition. Actually, even the data 
between these two regions are also scaled on a single 
curve 0 E 1 , and the exponent fl/v can be evaluated 
from this “nonequilibrium-to-equilibrium” scaling plot. 
Although this plot seems to have 3 parameters cr, j3/v and 
T m , fitting procedure is not so difficult. The data close to 
equilibrium almost only depend on ft/v, and there exists 
a constraint that the tangent of the short-time data is 
unity. Other exponents 7 /z/ and a/iz can be estimated 
from similar scaling plots of the magnetic susceptibility 
and specific heat. Then, the exponent v can be obtained 
from the hyperscaling relation, 2 — a = dv. 

However, direct evaluation of a/v from size depen¬ 
dence of the specific heat is not trivial in the 3D vector 
spin models, because the critical exponent a is known 
to be close to zero or even negative in these models. 
In such a case, the temperature derivative of magneti¬ 
zation El Ei, EH might be a good quantity. Since the 
spontaneous magnetization behaves as m s ~ (T c — T)& in 
the vicinity of T c , its temperature derivative behaves as 
dm s /dT ~ — (T c — T)^ _1 , and the size dependence of this 
quantity at T c is given by dm s (L;T c )/dT ~ — L( 1 “< 3 )/ 1 '. 

In order to evaluate the critical exponents precisely 
enough, the critical temperature T c should also be eval¬ 
uated precisely beforehand. In the traditional NER 
scheme based on the local-update algorithm, T c is evalu¬ 
ated from linearity of the log-log plot of relaxation data 
for a single system size in the power-law region. This 
region becomes fairly wide when rather large systems are 
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FIG. 1: (Color online) Nonequilibrium-to-equilibrium scaling 
plot of magnetization at T c = 1.442987J/A;b for L = 200, 280, 
400 and 560 with j3/v = 0.516, a = 0.47 and T m = 2.51. All 
the data seem to be located on a single curve within this scale. 
The dashed line stands for the tangent 1 as guides for eyes. 


considered. Usually, saturation to equilibrium behaviors 
is never observed in the standard NER analysis. On the 
other hand, the stretched-exponential region is rather 
narrow in the NER analysis based on the cluster algo¬ 
rithms. Especially in three or more dimensions, precise 
evaluation of T c from a single system size is challenged by 
saturation to equilibrium behaviors. Such attempts will 
be summarized in Appendix A fl2| . and the final-stage 
estimation of T c coupled with determination of critical 
exponents is explained in the present section. 


B. Evaluation of T c and fi/u from magnetization 

First, the scaling plot of magnetization © at the 
estimated critical temperature T c = 1.442987J//c B (to 
be explained later) with fj/v = 0.516, a = 0.47 
and T m = 2.51 is displayed in Fig. [lj where all the 
data points for L = 200 (averaged with 320, 000 ran¬ 
dom number sequences (RNS)), 280 (160,000 RNS), 400 
(80,000 RNS) and 560 (80,000 RNS) seems to be lo¬ 
cated on a single curve, as expected. Evaluation pro¬ 
cess of T c and p/u is shown in Fig. [2] where the data 
for T = 1.442984J/fc B , 1.442985 J/k B , 1.442987J/fc B , 
1.-442989 J/ks and 1.442990J/fc B (from bottom to top) 
in the saturation region are plotted with Eq. ®, and 
the data are almost on a single curve after the tun¬ 
ing of parameters. The finest scaling is observed at 
T = 1.442987J//c B with (3/v = 0.516, while deviations 
of the data at T = 1.442984J/fc B with P/v = 0.511 and 
at T = 1.442990J/fc B with p/v = 0.519 are not negli¬ 
gible. Since the direction of deviations is not system¬ 
atic as varying the system size, these deviations cannot 
be reduced anymore. We take similar scaling plots for 
1.442985J/fc B < T < 1.442989J/fc B , and observe com¬ 
parable scaling behaviors with gradually-varying expo¬ 



FIG. 2: (Color online) Enlarged plot in the vicin¬ 

ity of equilibrium corresponding to Fig. at T = 
1.442984J/fc B , 1.442985 J/k B , 1.442987 J/k B , 1.442989 J/k B 
and 1.442990J/&B with 0/v = 0.511, 0.512, 0.516, 0.518 and 
0.519, respectively. Error bars are smaller than symbols. 

nent 0.512 < P/v < 0.518. The fine scaling behaviors 
for T = 1.442985J/fc B < T < 1.442989J/fc B becomes as 
bad as that at T = 1.442984J/fc B by changing P/v with 
±0.002. Combining these results, we can safely conclude 

T c = 1.442987(2) J/k B , P/v = 0.515 ± 0.005. (5) 

When we fit the data for L = 560 at T c directly with 
Eq. ®, we have <r ps 0.49, which is consistent with that 
of the 3D XY model [15} and the exponent a « 1/2 
might be common in all the vector spin models char¬ 
acterized by the second-order phase transition. In this 
fitting the initial ~ 50 MCS data are used by minimizing 
the residue per data Q. With this exponent the data be¬ 
tween the stretched-exponential and equilibrium regions 
are not scaled so well as those plotted in Fig. [T] due to 
higher-order correction of simulation time, and the devi¬ 
ation from cr = l/2 may compensate such higher-order 
correction numerically. 

Evaluation process of a is given in Fig. [3] which ex¬ 
hibits the average mutual variance S 2 in the scaling plot 
like Fig. □ at T c = 1.442987 J/k B for a = 0.460 (tri¬ 
angles), 0.465 (crosses), 0.470 (circles), 0.475 (x-marks) 
and 0.480 (squares). After the optimization of param¬ 
eters P/v and r m with given a for L = 200, 280, 400 
and 560, variance of the data between different L at the 
same scaled simulation time t* = ( t/r m ) a — log L d l 2 ~^l v 
(between the law data for one size and the linearly- 
interpolated ones for other sizes) is averaged in the width 
At* = 0.2 (for example, from t* = —5.0 to —4.8). 

Apparently, the variance is very small at the onset of 
ordering (f* ps —5) and in the vicinity of equilibrium 
(f» ~ 2), and has two main peaks on the edge of the 
stretched-exponential ordering (—3 < t* < —1; region A) 
and between the stretched-exponential and equilibrium 
regions (—1 < t* < 1; region B). For a = 0.48, the 
variance is the smallest in the region A while the largest 
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FIG. 3: (Color online) Average mutual variance of the scal¬ 
ing plot in Fig. [l] for a = 0.460 (triangles), 0.465 (crosses), 
0.470 (circles), 0.475 (x-marks) and 0.480 (squares). Overall 
variance seems to be the smallest for a « 0.47. 


in the region B, and vice versa for a = 0.46. For a = 0.47 
the variance takes medium values in the both regions and 
these two values are almost the same, and therefore the 
most suitable within these three values. Behaviors for 
a = 0.475 and 0.465 just locate in between. Optimized 
parameters are f)/v = 0.515 and r m = 2.31 for a = 0.46, 
and f3/v = 0.517 and r m = 2.74 for er = 0.48. That is, 
the variance of fi/v is much smaller than that of a. 

Calculations for the variance of the next order of a are 
affected by the precision of optimization of parameters 
or the method of interpolation (linear or up to higher 
orders). Such efforts may not be productive, because 
error bars in /?/V are smaller than those from the scaling 
analysis in Fig. [2] (±0.002) for each temperature. Then, 
here we evaluate this exponent as a = 0.47(1) and use 
this estimate for analyses of other physical quantities. 


C. Evaluation of 7 / 1 / from magnetic susceptibility 

Although evaluation of T c based on the magnetic sus¬ 
ceptibility is also possible, here we use T c based on mag¬ 
netization and estimate 7 jv similarly to the scheme used 
in Fig. [2] First, the magnetic susceptibility at T c in the 
SW algorithm shows the stretched-exponential ordering, 

X(t, L ) = Jd ■ Sj) ~ exp [+{t/T x ) a ] . ( 6 ) 

This physical quantity has no size-dependent prefactor 
because duplicated random-walk growth of clusters (~ 
L~( d / 2 )x 2 -j cance u ec i with bulk response (~ L d ). Since 
equilibrium value of this quantity is scaled as % c ~ L 7 /", 
we have the following scaling form similar to Eq. 0 , 



FIG. 4: (Color online) Nonequilibrium-to-equilibrium scaling 
plot of magnetic susceptibility at T c = 1.442987J/A:b for L = 
200, 280, 400 and 560 with 7 / 1 / = 1.970, a = 0.47 and r x = 
0.584. All the data seem to be located on a single curve, and 
the dashed line stands for the tangent 1 as guides for eyes. 



FIG. 5: (Color online) Enlarged plot in the vicinity of 
equilibrium corresponding to Fig. [4] at T = 1.442985 J//cb , 
1.442987J/fc B and 1.442989J/fc B with 7 / 1 / = 1.977, 1.970 and 
1.967, respectively. Error bars are smaller than symbols. 


The scaling plot at T c = 1.442987J//cb with 7 jv = 1.970, 
a = 0.47 and t x = 0.584 is displayed in Fig. [2 where all 
the data points for L = 200, 280, 400 and 560 seems 
to be located on a single curve again. The data in the 
vicinity of equilibrium is shown in Fig. [5j together with 
the data for T = 1 .442985 J/ks (with ~/jv = 1.977) and 
1.442989 J//cb (with 7 /^ = 1.967). Discrepancy from 
such scaling behaviors becomes clear by changing 7 /^ 
with ±0.002. Combining these results, we conclude 

7 /V = 1.972 ±0.007. ( 8 ) 


log( X (f, L)L~^) ~ ( t/T x Y - log ± 7 /± (7) 
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FIG. 6: Time dependence of the reversed temperature deriva¬ 
tive of magnetization p(t) at T c = 1.442987J/fcB for L = 560. 


D. Evaluation of u from temperature derivative of 
magnetization and a from scaling relation 

As explained in the Overview part of this section, the 
quantity n(t) = —d(\m(t)\)/dT might be suitable for 
evaluation of v. Using the Hamiltonian of the model 
© as the energy variable, this quantity is expressed as 

n(L) = -±((\m\H)-(\m\)(H)). (9) 

This expression is nothing but the subtraction between 
two nearly-equal quantities. Although the error bar is 
expected to be large in equilibrium, early-time behaviors 
from the perfectly-disordered state might be different, 
and law data of time dependence of /x(£) for L = 560 
averaged with 80, 000 RNS is displayed in Fig. [6] 

The quantity ( \m\'H) monotonically decreases from 
zero as the simulation time increases, and it reaches 

-7 x 10 6 at t = 225 MCS for L = 560. That is, p(t) 

is obtained from the subtraction between two quantities 
of order of ~ 10 5 times larger. We may call the error bar 
of n(t) rather “small” even for such severe calculations. 
This figure also reveals that the time dependence of n(t) 
is quite nontrivial: It takes negative values in initial sev¬ 
eral steps, increases rapidly, arrives at a maximal value 
and gradually decreases, shows upturn again and finally 
saturates toward equilibrium. It is not possible to decide 
which region corresponds to the stretched-exponential 
time evolution a prior , and we have attempted some trials 
and investigated which assumption results in fine scaling 
behaviors in a wide parameter region. 

We finally find that the stretched-exponential decay 
before the upturn seems reasonable, 

p(t, L) ~ L d/2 exp[—(t/r^)' 7 ], (10) 

where the size dependence is derived by multiplying 
random-walk growth of clusters and bulk response. From 



FIG. 7: (Color online) Nonequilibrium-to-equilibrium scaling 
plot of the quantity p(t, L ) at T c = 1.442987J/fcB for L = 200, 
280, 400 and 560 with (1 — 0)/v = 0.889, a = 0.47 and = 
7.40. All the data except for the initial rapidly-increasing 
ones seem to be located on a single curve, and the dashed line 
stands for the tangent —1 as guides for eyes. 
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FIG. 8: (Color online) Enlarged plot in the vicinity of 
equilibrium corresponding to Fig. [7] at T = 1.442985 J/ks , 
1.442987 J/fcs and 1.442989 J/k B with (1-/3) /v = 0.890, 0.889 
and 0.888, respectively. Error bars are a few times larger than 
symbols, which exhibits large fluctuation of the data. 

Eq. (HOD and the size dependence in equilibrium, /x c (L) ~ 
we have the following scaling form, 

log ( f i(t,L)L-( 1 -W v ) ~ {t/r^Y - log l<*/ 2 -( 1 -0)/*'. 

( 11 ) 

The scaling plot at T c = 1. 442987 J/ks with (1 — (3)/v = 
0.889, a = 0.47 and = 7.40 is displayed in Fig. [7] 
where all the data points except for the initial rapidly- 
increasing ones for L = 200, 280, 400 and 560 seems to 
be located on a single curve. Actually, the stretched- 
exponential scaling region is rather narrow and precise 
evaluation of is difficult from these data. Then, we 
fit this parameter with the data after the upturn, and 
confirm the acceptable stretched-exponential scaling with 
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this estimate afterwards. The data in the vicinity of 
equilibrium is shown in Fig. [8] together with the data 
for T = 1.442985 J/k B (with (1 ~ 0)/v = 0.890) and 
1.442989J/fce (with (1 — f})/v = 0.888). Discrepancy 
from such scaling behaviors becomes clear by changing 
— with ±0.005, and we have (1 — 0)/v = 0.889(6). 
Combining this estimate with Eq. © , we arrive at 

l/v = 1.404 ±0.008 or v = 0.712 ± 0.004. (12) 

Furthermore, when this estimate is coupled with the hy¬ 
perscaling relation, 2 — a = dv, we have 

a = -0.136 ±0.012 or a/v = -0.192 ± 0.016. (13) 

E. Complementary evaluation of a/v from specific 
heat and v from scaling relation 

The specific heat is expressed as 

C(L) = «« 2 > - (H) 2 ) . (14) 

This expression is also the subtraction between two 
nearly-equal quantities, of order of ~ 10' times larger 
than C(L) itself in the vicinity of equilibrium for L = 560. 
Although this ratio is even larger than that of er¬ 

ror bars of C'(L) are smaller than those of n(L), because 
fluctuation of energy is much smaller than that of magne¬ 
tization. Although time dependence of this quantity is es¬ 
sentially the same as that of in Fig. [© the stretched- 
exponential decay of this quantity is quite apparent (as 
will be shown in Fig. [9]), 

C(t, L) ~ L d exp [—(t/rcY ], (15) 

where the size dependence originates from bulk response. 
Combining this scaling form with the equilibrium finite- 
size scaling C c (L) ~ L a l v , we arrive at 

log (C{t, L)L-“/ y ) ~ ( t/rcY - log L d ~ a /r (16) 

The scaling plot at T c = 1.442987J/fcB with a/v = 
0.074(5), cr = 0.47 and t c = 0.215 for L = 200, 280, 400 
and 560 is given in Fig. [9] The enlarged data after the 
upturn is shown in Fig.[lU]together with the data for T = 

1.442985J/Zcb and 1.442989J/fc B with a/v = 0.077(6) 
and 0.070(5), respectively. Error bars of the data are 
a bit larger than symbols, and error bars of estimates 
of the critical exponent a/v are also larger than those in 
other exponents. Combining these results we have a/v = 
0.074±0.009, and v = 0.651±0.002 from the hyperscaling 
relation. This estimate is not at all consistent with the 
one in Eq. m , and such discrepancy had already been 
discussed in the 3D Heisenberg model (l^ - l20j . 

In equilibrium simulations, such a puzzle of the spe¬ 
cific heat is known to be resolved by dividing it to the 
regular and scaling parts as C(L) = C reg — aL a / v fl6i ]. 



FIG. 9: (Color online) Nonequilibrium-to-equilibrium scaling 
plot of specific heat at T c = 1.442987J//cb for L = 200, 280, 
400 and 560 with a/v = 0.074, a = 0.47 and tc = 0.215. All 
the data except for initial several steps and those during the 
upturn seem to be located on a single curve, and the dashed 
line stands for the tangent —1 as guides for eyes. 



FIG. 10: (Color online) Enlarged plot in the vicinity of 
equilibrium corresponding to Fig. [5] at T = 1.442985 J/ks , 
1.442987 J/fcB and 1.442989 J/fc B with a/v = 0.077, 0..074 
and 0.070, respectively. Error bars are a bit larger than sym¬ 
bols. This figure exhibits fluctuation of the data near equilib¬ 
rium and large deviation of the data during the upturn. 

In the present case, we should also consider time depen¬ 
dence of C(L), namely the stretched-exponential scal¬ 
ing given in Eq. m- Now we interpret this form as 
an alternative definition of the scaled simulation time 
t* = ( t/rcY — log L d , and the parameter tc is deter¬ 
mined as the minimum of C(t,L) to be independent of 
system size as shown in Fig. Illf a'I. which corresponds to 
t c « 0.27. 

Now we assume that the regular part of the specific 
heat is independent of simulation time t. while the coef¬ 
ficient of the scaling part depends on t as 

C{t,L) = C leg -a(t)L a ^. 


(17) 
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FIG. 11: (Color online) (a) Specific heat at T c = 
1.442987J//cb for L = 200, 280, 400 and 560 plotted versus 
the scaled simulation time t* = ( t/rcY — log L d derived from 
the stretched-exponential scaling form m- The parameter 
tc ~ 0.27 is determined as the minimum of the specific heat 
becomes independent of L (emphasized by the dashed line), 
(b) Scaling of the specific heat including the regular part (11711 . 
Assuming a/v = —0.192 obtained from the scaling behavior 
of fi(t) in Eq. (11311 and tuning Greg ~ 5.05&B, the coefficient 
of the scaling part of the specific heat becomes independent 
of L in the whole region satisfying C re g > C(i, L). 


The point of this formulation is that C reg > C(t,L) or 
a(f) > 0 is satisfied in the scaling region, and that the 
exponent a can be negative even though C{t, L) increases 
as the linear size L increases. Although it is possible to 
estimate C reg and a/v from the scaling plot based on 
Eq. (TT71) in principle, error bars of the estimates must be 
much larger than that of Eq. m because there exists 
an extra fitting parameter C reg while the fluctuation of 
the data is comparable to that of /). 

Then, we only evaluate Creg by assuming the exponent 
given in Eq. (1131) and confirm the validity of the scaling 
(fTTll . In Fig. HHb), the quantity (Creg — C(t, V))L~ a l v 
is plotted versus the scaled simulation time with C reg ~ 
5.05/cb at T c = 1.442987 J/k B for L = 200, 280, 400 and 
560, namely the ^-dependence of the coefficient of the 


Refs. 


Method 

P 

7 

V 

present I 


MC 

0.367(4) 

1.404(9) 

0.712(4) 

present II 


MC 

0.367(2) 

1.403(6) 

0.712(3) 

[16] (1993) 


MC 

0.362(3) 

1.389(14) 

0.704(6) 

[21] (2001) 


MC 

0.3685(11) 1.393(4) 

0.710(2) 

[22] (2002) 


MC 

0.3691(6) 

1.3957(22) 0.7113(11) 

[22J (2002) 

MC+HTE 

0.3689(3) 

1.3960(9) 

0.7112(5) 

[23] (1997) 

HTE (sc) 

0.3710(13) 1.406(3) 

0.716(2) 

[23] (1997) 

HTE (bee) 

0.3700(13) 1.402(3) 

0.714(2) 

[24] (1998) 


e-Exp. 

0.367 

1.39 

0.708 

[25] (1998) 


e-Exp. 

0.3655(35) 1.382(9) 

0.7045(55) 

[25] (1998) 

d 

= 3 Exp. 

0.3662(25) 1.3895(50) 0.7073(35) 

[26] (2001) 

d 

= 3 Exp. 

0.3655(5) 

1.3882(10) 0.7062(7) 


TABLE I: List of critical exponents estimated in the present 
study and in several previous numerical studies. 


scaling part a(t). Apparently, the scaling behavior in 
Fig. [TU b) seems much better than that in Fig. [10] or 
the improved scaling form mu with negative a might be 
better than the naive scaling form (11611 with positive a. 


IV. DISCUSSION 

First, combining the scaling relation a + 2(3 + 7 = 2 
and the hyperscaling relation, we have 2/3/iz + 7 / 1 / = d. 
On the other hand, from Eqs. © and © we result in 
2/3/V+ 7 /;/ = 3.002±0.012. That is, the critical exponent 
a/v can be evaluated from the temperature derivative 
of the magnetization without the hyperscaling relation, 
which is not so trivial as other scaling relations. 

Next, we compare the estimates of /3, 7 and is obtained 
from Eqs. ©, © and (fl^l) (present I) with several previ¬ 
ous studies based on the Monte Carlo method 0 [piEl, 
high-temperature expansion [22I [23j , e-expansion |24l.l25l| 
or d = 3 expansion [25], |26j] in Table [Q More comprehen¬ 
sive list of references was given in Ref. [22|. Note that 
our estimates shown above are based on “safe” evalua¬ 
tion. That is, the critical exponents are determined in 
order to cover the whole possible range of the critical 
scaling. If we assume that T c locates in the middle of 
the scaling region and the critical exponents are evalu¬ 
ated just at T c = 1.442987(1) J//cb, reasonable reduction 
of error bars is possible (present II). 

From this table we may at least conclude that the esti¬ 
mates of critical exponents in the present study are con¬ 
sistent with those in previous studies, even though er¬ 
ror bars are not so small. The dominant origin of sta¬ 
tistical errors in our estimates is large fluctuations of 
li(t,L ). Since all the critical exponents evaluated with 
finite-size scaling analyses are scaled by is, Precision of 
this exponent is crucial. Nevertheless, the main purpose 
of the present study is to confirm the nontrivial stretched- 
exponential critical scaling in the 3D Heisenberg model, 
and evaluation of the critical exponents is not more than 
consistency check of this alternative scheme. 

Actually, information of the critical exponents is not 


















included in the stretched-exponential critical scaling of 
physical quantities in itself. Such information is intro¬ 
duced from the finite-size scaling of physical quantities in 
equilibrium, and the nonequilibrium-to-equilibrium scal¬ 
ing scheme enables to extract such information from the 
off-scaling data prior to equilibrium behaviors. The most 
promising application of the present scheme is to charac¬ 
terize phase transitions with nonequilibrium behaviors. 
In the standard NER scheme based on the local-update 
algorithms, power-law behaviors take place in general 
even in the BKT j27| or first-order [28| phase transitions, 
and such characterization was made in a different manner 
by assuming the nature of phase transitions in advance. 
The situation might be different in the present cluster 
NER scheme, and investigations along this direction is 
now in progress il5[ . 

The stretched-exponential nonequilibrium critical scal¬ 
ing initially found in the 2D Ising model Q is now con¬ 
formed in the 3D Heisenberg model. The exponent of 
stretched-exponential time dependence a ~ 1/2 seems 
consistent with that of the 3D XY model |U|, while 
not consistent with that of the 2D or 3D Ising models, 
a « 1/3 @, 0 . Even though our simulations on 3D vec¬ 
tor spin models are based on the “embedded-Ising-spin” 
formalism, the exponent a is different from that of the 
Ising models. That is, this exponent is not determined 
by mere formalism, but by fundamental physical prop¬ 
erties. Further details of such universal behaviors of the 
exponent a will be explained elsewhere fl3[ . 


V. SUMMARY 

Critical behaviors of the 3D Heisenberg model are in¬ 
vestigated with the nonequilibrium ordering from the 
perfectly-disordered state in the Swendsen-Wang (SW) 
algorithm characterized by the stretched-exponential 
critical scaling form. Calculations from the perfectly- 
ordered state in the SW algorithm result in larger initial¬ 
time deviation and slower off-critical deviation, and those 
from the perfectly-disordered state in the Wolff algorithm 
show a power-law behavior similar to that in local-update 
algorithms. From simulations based on the “embedded- 
Ising spin” formalism we have the nonequilibrium-to- 
equilibrium scaling similarly to the 3D XY model, and 
all the data for L = 200, 280, 400 and 560 seem to be on 
a single curve in the whole simulation-time region with 
the stretched-exponential exponent a = 0.47(1). 

Precise values of the critical temperature and the crit¬ 
ical exponent /v are evaluated simultaneously from the 
scaling analysis of the magnetization. Similarly, the crit¬ 
ical exponent 7 / 1 / is obtained from the scaling analysis 
of the magnetic susceptibility. The critical exponent v 
is evaluated from the scaling analysis of the tempera¬ 
ture derivative of magnetization, and the critical expo¬ 
nent a < 0 can be derived from the hyperscaling relation. 
A consistent estimate of ajv can be obtained from the 
specific heat by assuming the scaling form with the size- 


independent regular part and the negative scaling part. 

Although information on standard critical behaviors 
such as the critical exponents is not included in the 
stretched-exponential critical scaling behaviors, such in¬ 
formation can be extracted from off-scaling behaviors 
with the nonequilibrium-to-equilibrium scaling scheme, 
and the critical exponents comparable to previous stud¬ 
ies can be obtained from the present framework without 
full equilibration. 
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Appendix A: Evaluation of critical temperature 
from single system 

1. Double-log plot 

The simplest way to evaluate the critical temperature 
T c from the early-time dependence of magnetization © 
or © is the double-log plot. When both sides of these 
equations are divided by the coefficient of rhs (expressed 
as C) and take logarithm twice, we have 

log(log((|m(t)|)/C')) - =F0-log(i/T m ). (Al) 

Since Eq. © is obtained by relaxation from the perfectly- 
ordered state, the coefficient C is of order of unity, and a 
naive plot assuming C = 1 may be satisfactory enough. 
On the other hand, Eq. © is obtained by ordering from 
the perfectly-disordered state, and the small coefficient 
C should be taken into account explicitly. 

Practically, the unknown parameter C can be deter¬ 
mined by aligning the early-time data on a straight line, 
as shown in Fig. fl2l' a') with the data for L = 100 with 
2,000 RNS in the SW algorithm from the perfectly- 
disordered state. Apparently, the ordering process in 
the double-log plot is faster or slower than linear for 
T < lAAJ/ks or T > 1.45J/fcB, respectively, which 
means 1.44J//cb < T c < 1.45J/fcB- Tangent of the 
dashed line is 0.5, which is consistent with a « 0.5. Then, 
the variance of temperatures is reduced of one order in 
Fig. []©b) with the data for L = 200 with 2, 000 RNS, 
which results in T c ss 1.443J/fcs with cr ~ 0.5. 

Differently from the standard NER scheme based on 
the power-law behavior of physical quantities, the present 
scheme requires the rescaling of physical quantities with 
the unknown coefficient of stretched-exponential behav¬ 
iors, and therefore precise evaluation of T c is difficult. 
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t [MCS] 



t [MCS] 


FIG. 12: (Color online) Double-log plot of magnetization in 
the SW algorithm from the perfectly-disordered state for var¬ 
ious temperatures around T c for (a) L = 100 and (b) L = 200. 
The dashed line corresponds to the stretched-exponential or¬ 
dering with <7 = 0.5. 
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2. Semi-log plot with fixed a 

If we intend to estimate more accurate T c from the 
double-log plot, we may determine the linearity of the 
early-time data assuming the exponent a = 0.5. How¬ 
ever, as long as er is fixed a prior , plotting log(|m(t)|) 
versus t a is more straightforward and precise, because 
no adjusting parameter other than cr is included. Such 
semi-log plot with cr = 0.5 is shown in Fig. [TUI a) for 
L = 560 with 10,000 RNS around T = 1.443 J/k B . Eval¬ 
uation of T c is not an easy task because the data start 
deviating at t = 50 ~ 60 MCS, while they already be¬ 
gin to saturate at t = 70 ~ 80 MCS. Observing the 
data on the onset of deviation carefully, we may conclude 

1.4429J/Zcb < T c < 1.4430J/fes- 

Then, the variance of temperatures is reduced of one 
more order in Fig. 0b) with the data for L = 800 
with 5, 000 RNS. Although the data seem to suggest 
1.44294 J//cb < T c < 1.44295J/fcB, this estimate of T c is 


FIG. 13: (Color online) Semi-log plot of magnetization versus 
t" with o = 0.5 in the SW algorithm from the perfectly- 
disordered state for various temperatures around T c for (a) 
L = 560 and (b) L — 800. The dashed line corresponds to T c 
obtained from these figures. Discrepancy with the estimate 
in Section III reveals the limitation of this graphical scheme. 


not consistent with the scaling plot given in Eq. (U) at all. 
Even though the data starts deviating at t = 80 ~ 100 
MCS, they already begin to saturate in such simulation¬ 
time region. Downward-bending behavior due to satura¬ 
tion and upward-bending behavior below T c are cancelled 
with each other in this region, and therefore T c is under¬ 
estimated when it is determined by the linearity of the 
nonequilibrium data. The scaling analysis explained in 
Section III is more reasonable and precise. 
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FIG. 14: (Color online) Semi-log plot of magnetization versus 
t a with a = 0.5 in the SW algorithm from the perfectly- 
ordered state for various temperatures around T c for L = 
560. The dashed line emphasizes the region in which the 
stretched-exponential decay holds well, and it clarifies early- 
time rounding and saturation of data earlier than deviation 
induced by varying temperatures. 


Appendix B: Nonequilibrium relaxation from the 
perfectly-ordered state 

1. Numerical results in the SW algorithm 

Relaxation data of magnetization from the perfectly- 
ordered state in the SW algorithm for L = 560 with 400 
RNS are displayed in Fig. [lj] with a semi-log plot versus 
t a , a = 0.5 similarly to Fig. [13] This figure shows mer¬ 
its and demerits of this scheme. First, although average 
number of RNS is one order smaller than that in Fig. 
[T3l diversity of the data looks comparable. On the other 
hand, stable stretched-exponential decay only starts from 
t ~ 20 MCS as emphasized by the dashed line, while 
such behavior is observed from the beginning in Fig. 1121 
Moreover, the onset of deviation of the data induced by 
variance of temperature is t ~ 80 MCS, which is later 
than that in Fig. [T3] and even later than the onset of 
saturation in Fig. [TT] f ~ 70 MCS, which is compara¬ 
ble to that in Fig. [13] That is, evaluation of T c from 
deviation of the data does not make sense, because this 
deviation is already affected by saturation behavior. 

2. Numerical results in the Wolff algorithm 

Relaxation data of magnetization from the perfectly- 
ordered state in the Wolff algorithm for L = 560 with 
100 RNS are shown in Fig. [15] with a standard semi-log 
plot versus t, or the critical magnetization decays ex¬ 
ponentially in this case. Similarly to the corresponding 
relaxation in the SW algorithm, stable exponential de¬ 
cay only starts from t ~ 20 MCS, and saturation of the 



FIG. 15: (Color online) Semi-log plot of magnetization in the 
Wolff algorithm from the perfectly-ordered state for various 
temperatures around T c for L = 560. The dashed line empha¬ 
sizes the region in which the simple exponential decay holds 
well. Although time dependence looks similar to that in Fig. 
ELD we take a = 1 in this figure. 


data begins earlier than deviation induced by variance 
of temperature. Rather small number of RNS is due to 
larger numerical cost than that in the SW algorithm. Al¬ 
though required simulation time for equilibration looks 
much smaller than that in the SW algorithm, this time 
is scaled by the flipped cluster size in each step in the 
Wolff algorithm. Since the exponent /?/V is rather large 
(s=s 0.5) in the present model, the critical magnetization 
m c (L) rapidly decreases as the linear size L increases. 
This situation is quite different from that in the 2D Ising 
model Q, where the small exponent f3/v = 1/8 ensured 
efficient update in the Wolff algorithm. 

Aside from numerical efficiency, it is theoretically in¬ 
teresting that the exponent of the stretched-exponential 
decay (simple exponential decay can be regarded as the 
ct = 1 case) depends on the species of the cluster algo¬ 
rithm. Such behaviors did not take place in the 2D Ising 
model Q, where the decay exponent from the perfectly- 
ordered state is the same in the SW or Wolff algo¬ 
rithms. (Note that ordering behaviors from the perfectly- 
disordered state is quite different in the two algorithms 
even in the 2D Ising model [§,].) Although the direction 
of flip is trivial in the Ising model, it is chosen randomly 
in the Heisenberg model. In each step the alignment is 
common in the whole system in the SW algorithm, but it 
is random from cluster to cluster in the Wolff algorithm. 
Absence of correlation between clusters is characteristic 
to behaviors above T c , and it results in the exponential 
decay. That is, the critical slowing down does not exist 
in the relaxation from the perfectly-ordered state in the 
Heisenberg model in the Wolff algorithm. Although this 
behavior is favorable for equilibration, it is not suitable 
for NER analyses. 
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